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ABSTRACT 

Disintegration of sunspots (and starspots) by fluxtube erosion, originally pro¬ 
posed by Simon & Leighton, is considered. A moving boundary problem is formu¬ 
lated for a nonlinear diffusion equation that describes the sunspot magnetic held 
prohle. Explicit expressions for the sunspot decay rate and lifetime by turbulent 
erosion are derived analytically and verihed numerically. A parabolic decay law 
for the sunspot area is obtained. For moderate sunspot magnetic held strengths, 
the predicted decay rate agrees with the results obtained by Petrovay & Moreno- 
Insertis. The new analytical and numerical solutions signihcantly improve the 
quantitative description of sunspot and starspot decay by turbulent erosion. 

Subject headings: dihusion — turbulence — Sun: magnetic helds — sunspots — 
stars: magnetic held — starspots 


1. Introduction 

Bumba (1963) investigated how the areas of large, slowly decaying sunspots decrease 
with time. His data analysis suggested that the sunspot area A decreases linearly with time 
t: 

A{f) = Ao- At, (1) 

where the decay rate A is a constant. The result is consistent with the Gnevyshev-Waldmeier 
relation T ~ Aq, where T is the sunspot lifetime and Aq is its initial area (see, e.g., Petrovay & 
van Driel-Gesztelyi (1997) for a review). Following Bumba (1963), sunspot observations were 
usually interpreted in terms of the linear decay law for the sunspot area (e.g., Robinson & 
Boice 1982). Yet it is difficult to distinguish linear and nonlinear decays observationally, and 
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observations have also been interpreted using a parabolic decay law, with A{t) a decreasing 
quadratic function of time (Moreno-Insertis & Vazquez 1988; Martinez Fillet et ah 1993). 

On the theoretical side, Meyer et ah (1974) argued that the linear decay law is a con¬ 
sequence of turbulent diffusion of the magnetic field across the whole area of a sunspot and 
expressed the constant decay rate in terms of a constant uniform diffusivity (see also Krause 
& Riidiger 1975). 

Simon & Leighton (1964) inferred from observations that the gradual disintegration of 
sunspots is due to “erosion” of the penumbral boundaries by supergranular flows, which 
occurs when bits of magnetic held are sliced away from the edges of the sunspot and swept 
to the supergranular cell boundaries. In contrast to the model of Meyer et ah (1974), such 
erosion can occur if the turbulent diffusivity associated with the hows is suppressed within the 
spot (Petrovay & Moreno-Insertis 1997). Alternative theoretical approaches were reviewed 
by Solanki (2003). 

Petrovay & Moreno-Insertis (1997) developed the turbulent erosion model mathemati¬ 
cally, taking into account the dependence of the turbulent dihusivity on the magnetic held 
strength. The dihusivity rapidly decreases if the magnetic held exceeds an energy equiparti- 
tion value (Kitchatinov et al. 1994). As a result, a current sheet is formed around the spot. 
The model leads to the parabolic decay law, specihed by a constant inward speed w of the 
current sheet, viz., 

A{t) = 7r(ro - wtf (2) 

for a circular hux tube (sunspot) of an initial area Aq = Moreover, the model yields w ~ 
1/ro, and so it agrees with the Gnevyshev-Waldmeier relation. Petrovay & Moreno-Insertis 
(1997) concluded that solar observations are consistent with turbulent erosion based on a 
granule-size dihusion length. Petrovay & van Driel-Gesztelyi (1997) presented observational 
evidence in favor of the parabolic decay rate, predicted by the turbulent erosion model, 
although an independent magnetohydrodynamic simulation suggested that the sunspot decay 
law is almost linear (Riidiger & Kitchatinov 2000). Petrovay et al. (1999) also explored the 
ehect of a preexisting “plage” held on the decay rate, whereas Ghatterjee et al. (2006) applied 
the model to the development of twist in a hux tube rising through the solar convection zone. 

The analytical results of the turbulent erosion model have been recently used to comple¬ 
ment numerical simulations of sunspot formation and decay (e.g., Rempel & Gheung 2014). 
The model has also been applied to starspots, with a goal of using the starspot decay data to 
place constraints on the magnetic dihusivity, which may be useful for dynamo models (e.g., 
Strassmeier 2009; Bradshaw & Hartigan 2014). 

It is worthwhile to revisit the turbulent erosion model of sunspot decay. The original 
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calculation of Petrovay & Moreno-Insertis (1997) was guided by numerical results and one¬ 
dimensional analytical solutions. A dimensional argument was used to estimate the magnetic 
field gradient at the sunspot edge: 


( 3 ) 


dB Be 

dr tq ’ 

where Bg is the magnetic held value above which the turbulent diffusivity is assumed to 
be suppressed (see equations (11) through (16) in Petrovay & Moreno-Insert is 1997). In 
addition, their numerical estimate for the sunspot lifetime appears to be based on an estimate 
of the current sheet speed w rather than on direct computation. 


A rigorous derivation of the sunspot decay law is necessary if the theory is to be used 
to develop reliable predictive tools. Explicit analytical predictions of the turbulent erosion 
model could complement more detailed numerical (e.g., Hurlburt & DeRosa 2008) and em¬ 
pirical (Gafeira et al. 2014) models of sunspot decay. Hence our aim is to put the turbulent 
erosion model on a hrmer footing. We do this by formulating a moving boundary problem 
(Carslaw & Jaeger 1959; Crank 1984) for the model and solving it to derive a prediction 
for the sunspot decay law. In the remainder of the paper, we present the new analytical 
(Section 2) and numerical (Section 3) results and their discussion (Section 4). 


2. Formulation of the problem and analytical results 


In order to model the turbulent erosion of a sunspot, we follow Petrovay & Moreno- 
Insertis (1997) and consider the evolution of a cylindrically symmetric magnetic flux tube. 
The magnetic held B = B{r,t)z is described by the diffusion equation 


^ (rD^] . 

dt r dr \ dr J 


( 4 ) 


Here t is time, and r is the distance from the ; 2 -axis. 

The turbulent diffusivity D = D{B) is strongly suppressed when the magnetic held 
exceeds an energy equipartition value Be = y/Anpu where p is the mass density and u 
is a characteristic turbulent speed (e.g., Kitchatinov et al. 1994). For instance, taking a 
photospheric value of p ~ 2 x 10“^ g cm“^ and a granular value of u ~ 2 x 10^ cm s“^ yields 
Be ~ 400 G (Petrovay & Moreno-Insert is 1997). To simplify the analytical treatment, we 
assume 

D{B) = Dq = const, B < Be (5) 


and 


D{B) = 0 , B > Be. 


( 6 ) 
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The initial value problem is specified by the field profile 


B{r, 0) = -Bo = const, 0 < r < ro, 


(7) 


and .B(r, 0) = 0 otherwise (see Tlatov & Pevtsov (2014) for recent data on sunspot magnetic 
fields). We nondimensionalize the problem by measuring the magnetic fields, times, and 
distances in units of B^, tI/Dq, and ro, respectively. 

The sunspot size decreases with time because the magnetic flux is removed by diffusion. 
The strongly nonlinear dependence of the diffusivity D on the magnetic field strength leads 
to the formation of a tangential discontinuity at the edge r = re of the flux tube. Physically, 
the magnetic field discontinuity at re(t) corresponds to a current sheet at the sunspot edge, 
where the magnetic flux removal is made possible by a strongly localized electric current. 

It is useful to observe that the problem at hand is mathematically similar to a moving 
boundary problem in the theory of heat conduction, and so we can use existing methods of 
analysis. In particular, an analog of the Stefan condition is obtained by the integration of 
the governing diffusion equation across the moving boundary r = re{t) (Carslaw & Jaeger 
1959). Allowing for the tangential discontinuity at r = re(t), we substitute 

B{r,t) = Bo + {B- Bo)H[r - r,(f)], (8) 


where H is the Heaviside step function, into equation (jl]) and integrate across the disconti¬ 
nuity (from Te — 0 to Te + 0). The result is 


{Bo- 



dB 


dr 


5 

r=re+0 


where we used Bije — 0, f) = .Bq, Bije -|- 0, t) = 1, and D{re — 0) = 0. 


(9) 


To find an approximate analytical solution, we use the pseudo-steady-state approxi¬ 
mation that can be adopted when the rate of change A is small compared with a global 
diffusion rate ~ 1, making it possible to neglect the term dB/dt in Equation (jl]). Physically, 
the magnetic field profile near a moving boundary relaxes to a pseudo-steady state on a time 
scale 6tE> — {SrY/D where 5r ~ -feht is the displacement of the boundary re(t) in a time 
6t. The approximation is valid if the relaxation is sufficiently rapid, say ii SId -C 5t. In 
our dimensionless variables, we have Jr ~ rg < 1 and .0 = 1, and so ~ -fe. If T is 

the sunspot lifetime, we use rg ~ T~^ to infer that, as long as T 1, the approximation is 
globally valid in the range 1 < f < T — 1. We show below that roughly T ~ .Bq — 1. Con¬ 
sequently, the pseudo-steady-state approximation becomes more accurate as Bq increases. 
Detailed analysis of the accuracy of the approximation can be found in standard textbooks 
on heat conduction (e.g.. Crank 1984; Hill & Dewynne 1987). 
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Inside the spot, the vanishing diffnsivity implies that the magnetic held is constant: 


B{r < re(t),t) = Bq. 

(10) 

Ontside the spot, the psendo-steady-state held satishes 


l_d_ ^ Q 

r dr \ dr J ’ 

(11) 

and so 

dB 

r—— = const. 
dr 

(12) 

Eqnation ((5|) gives the bonndary condition 


B{r = re{t),t) = 1, 

(13) 


which wonld he B = B^, in dimensional nnits. The magnetic held dihnsion ontside the spot 
canses the held to become negligibly small at some r = rf{t) ontside the spot. Solntions of 
the standard dihnsion eqnation in two dimensions snggest that rf{t) = (Carslaw & 

Jaeger 1959). Thns we set 

B{r = = 0. (14) 

The solntion of eqnation ([12]), satisfying the bonndary conditions at rg and Vf, is given by 


B(r > re{t),t) 

On snbstitnting this into eqnation (|9|), we get 


ln(rV2f) 
ln(r2/2t) ■ 


4 


(15) 


(16) 


We also obtained a similar diherential eqnation for re(f) nsing an independent heat- 
balance approximation (e.g.. Crank 1984). We do not present the resnlts here: althongh 
the approach reqnires longer calcnlations, it does not appear to be more accnrate than the 
psendo-steady-state approximation. 

Eqnation (ITHil does not appear to have a solntion in elementary fnnctions. The mag- 
nitnde of its right-hand side is of order nnity, which yields an order-of-magnitnde estimate 
rl{t) ~ 1 — t/(i?o — !)• Conseqnently, we have T ~ i?o — 1 and r‘l{T/2) ~ 1/2. Next 
we obtain a more accnrate solntion of eqnation flTbD . An approximate polynomial solntion 
wonld be convenient for comparison with the available observational resnlts and theoretical 
predictions. We nse a qnadratic approximation: 

~ Co -|- Ciit — T/2) -f C2{t — T/2)^, 


(17) 




where 


( 18 ) 



i=r /2 


and 


1 


C2 = 


2 df^ 


t=Tl2 


(19) 


We expand rlit) abont t = T/2 because this is where the pseudo-steady-state approximation 
is expected to be most accurate. The remaining constants Cq and T are dehned by the 
conditions 



r2(0) = 1 

(20) 

and 


rl(T) = 0. 

(21) 

Equations flTT)). fl2Ul). and (1^ give 

rl{t) 

~ 1 -|- (ci — C2T)t C2t^■ 

(22) 

Here 


T = - 

(23) 

Cl 

is the sunspot lifetime, unless T' <T where 


rp! _ Cl 

C2 

(24) 


is the other root of the equation = 0. 


We evaluate the constant Ci by substituting the order-of-magnitude estimates t = T/2 ~ 
(To — l)/2 and rl{T/2) ~ 1/2 into equation flTHl) . This yields an accurate expression for ci 
because T/2 and r‘l{T/2) only appear in the argument of the logarithm in equation fllbp . 
The resulting prediction for the sunspot lifetime is as follows: 

T = i(To-l)ln2(To-l), (25) 

which should be compared with equation (16) in Petrovay & Moreno-Insertis (1997) for the 
inward speed w = —f-e of the current sheet. In their model, w = const and the sunspot 
lifetime is given by Tpm = Te/w, which leads to 


Tpm = 


(26) 


in our dimensionless variables. The same result (up to a numerical coefficient) is obtained 
by nondimensionalizing our equation ([3]), substituting it into equation ([9]), and assuming 
he = const. 
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Differentiation of eqnation ffT6|) with respect to time yields 


C2 = 


{Bo - l)[Hrll2t)Y 


t {Bo - l)rl\n{r‘^/2t) 


(27) 


t=T/2 


Again nsing t = T/2 ~ {Bo — l)/2 and r‘l{T/2) ~ 1/2 is jnstified when these quantities 
appear in the argument of the logarithm. Therefore, 


C2 = 


1 + 


ln2(Eo-l) L 2r2(T/2) 


1 

2^2 ’ 


(28) 


where T is dehned by equation fl2S]) . The solution below can be used to verify that (T/2) = 
l/2 + 0(l/ln2(i?o —1)). Thus using rg(T/2) ~ 1/2 in equation fl28|) only leads to a relatively 
small error of order l/[Tln2(i?o — 1)]^, and we get 


2 

T2 ln2(So - !)■ 


(29) 


Collecting the results, we obtain a parabolic decay law for the sunspot area: 

2/ N . A 2 \ t 2 

~ ^ In 2{Bo - 1) j T ^ In 2(5o - 1) 7^' 

The sunspot lifetime is given by T in equation fl2^ if Bo > B^ and by 



r=l(B„-l)lln 2 (B„-l)|^ 

(31) 

if 1 < i?o < S*, where 

corresponds to T = T'. 

B, = l + eV2 ^ 4.7 

(32) 


Our explicit analytical solution for r^{t) provides an improved quantitative description 
of sunspot decay by turbulent erosion. Notably, if Bo = B^,, our solution predicts a constant 
decrease rate w = 2/{B^ — 1) 0.54 of the fluxtube radius: 


re{t) 1 — wt, 


(33) 


as in the parabolic decay law, predicted by Petrovay & Moreno-Insertis (1997). More gener¬ 
ally, we obtain a constant speed approximation 


(2 + lii2(B„-l)) T 


( 34 ) 















by defining w = —re(0) in our solution. If A{t) = Tir^ is the sunspot area, the accuracy of 
the approximation can be quantified by calculating 

2A 

which would be unity in the model 

On returning to the original dimensional quantities, we get the lifetime-size scaling 
T ~ Aq, where Aq = vrrg is the initial cross-sectional area of the flux tube. This result 
formally agrees with the Gnevyshev-Waldmeier relation for sunspot lifetimes. It is worth 
stressing that the statistical nature of the relation should follow from the strong dependence 
of T on the spot magnetic field Bq. 


t=o 


81n2(Eo - 1) 

|2 + lii2(B„-ip' 


(35) 


of Petrovay & Moreno-Insertis (1997). 


3. Numerical results 

The analytical results obtained in Section [2] may be tested by numerical solution of 
equation (0]). Following Petrovay & Moreno-Insertis (1997), we assume the analytical forms 
for the diffusivity and the initial field profile: 

= ( 37 ) 

where we use the non-dimensionalisation introduced in Section |2l The parameter in 
equation fl5^ determines the strength of the suppression of diffusion by the field, and the 
parameter in equation (137)1 specifies the initial spot profile. In the following we choose 
cxb = 22 to model an isolated flux tube with nearly constant internal field strength, and 
an = 7, to represent strong suppression of diffusion. For the purpose of numerical solution, 
the radius of the spot at time t is defined by the condition 

S(re,f) = ^5o. (38) 

We solve equation (|1]) using a Crank-Nicolson scheme (e.g. Press et ah 1992) which is 
described in the Appendix. The diffusion equation is evolved in time in the region 0 < r < 
with the boundary condition dB /Sr = 0 at r = 0 and with a boundary condition at r = 
which allows loss of flux from the region. Note that Petrovay & Moreno-Insertis (1997) used 
a less realistic condition dB/dr = 0 at an outer boundary (at r = 10), and their numerical 
solution was based on a Lax-Wendroff scheme. 
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Figure [T] illustrates the numerical solution for the case Bq = 7. The solid curves in the 
hgure show the numerical result for B{r, t) as a function of r for times t = 0, t = 0.5T, and 
t = 0.95T, where T is the analytical decay time, dehned by equation fl2^ . The solutions are 
shown for the region r < where = 7 is the outer boundary of the numerical domain. 
Figure [T] also shows the analytical solution at the same times, following equations ffTOj) 
and flT^ with the spot radius dehned by equation fl30|) . The spot decays more rapidly in 
the analytical solution, and the magnetic held outside the spot decreases more rapidly with 
increasing radius. The numerical solutions illustrate how the initial central hux concentration 
is redistributed to larger radius by dihusion, leading to an initial increase in held strength at 
points external to the spot. The qualitative behavior of the numerical solution is generally 
well reproduced by the analytical solution. 



Fig. 1.— Magnetic held versus radius at times f = 0, f = 0.5T, and t = 0.95T, for the case 
Bq = 7. The solid curves show the numerical solutions (for the parameter choices an = 7, 
as = 22), and the dashed curves show the analytical solutions. 

Figure [2] shows the square of the sunspot radius as a function of time for the same case 
Bq = 7. The solid curve shows the numerical solution, with Ve dehned by equation fl38l) . 
and the dashed curve shows the analytical solution dehned by equation fl30|l . The analytical 
solution decays more rapidly than the numerical solution, but both clearly show the departure 
from a linear decay law. The analytical estimate for the decay time is T = 3.72, and the 
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numerical decay time is 5.02. 

Figure E] plots numerically determined sunspot decay times versus central field strength 
Bq (crosses). The decay time is seen to depend almost linearly on field strength. The dashed 
curve shows the analytical results of Section |2l Recall that equation fl30|) defines two times at 
which re{t) = 0, namely T in equation (|25|) and T' in equation (El]). The decay time for the 
spot is given by T if Bq > and by T' if Bq < B^, where B^ is defined by equation fl3^ . 
The dotted vertical line in Figure El indicates the threshold value R*. Figure El also shows the 
decay time in the Petrovay & Moreno-Insertis (1997) model. Our analytical predictions agree 
with the numerical results: in particular the rates of increase of decay time with Bq are quite 
similar. Although our analytical model underestimates the decay times, it is significantly 
more accurate than the earlier model. 

Finally, we emphasize that our calculation generally yields a time-dependent rate of 
decrease of the fiuxtube radius re(t). Equation fl30|) predicts that the deviation from the 
parabolic decay law, derived by Petrovay & Moreno-Insertis (1997), should increase as the 
initial magnetic field Bq increases. As a result, a linear decay law (rather than a parabolic 
one) should become more accurate as Bq increases, although the logarithmic dependence on 
Bq makes the effect rather weak. Figure E] shows the effect of doubling the field strength 
Bq on the shape of the function rl{t). While the computation time and numerical errors 
increase for larger Bq, we do see numerical evidence that the decay law becomes more linear 
for a larger initial magnetic field, which is consistent with our analytical prediction. 


4. Discussion 

We have presented in this paper a quantitative theory of sunspot decay by turbulent 
erosion, considered as a moving boundary problem. The physical mechanism of sunspot 
erosion was proposed by Simon & Leighton (1964), and a sunspot decay law due to turbulent 
erosion was derived by Petrovay & Moreno-Insertis (1997) (see also Petrovay et al. 1999 
and references therein). Although Petrovay and collaborators correctly identified the key 
dependence of the decay rate on the sunspot magnetic field Bq, the accuracy of the analytical 
predictions was limited: for instance, we have shown that the numerically computed sunspot 
lifetime is about a half of that predicted. 

Our equation fl25p for the sunspot decay time T is an improvement on equation fl26p . 
derived by Petrovay & Moreno-Insertis (1997). Equation fl30P confirms that the decay law 
for the sunspot area A{t) = vrUg is in general parabolic, as long as higher-order terms in t/T 
can be neglected. Equation fl5Up also quantifies the accuracy of the assumption, made by 
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Fig. 2.— Sunspot radius squared versus time, for the case shown in Figure 1. The solid 
curve shows the numerical solution, and the dashed curve shows the analytical solution. 
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Fig. 3.— Decay time versus sunspot field strength Bq. The crosses indicate results for 
numerical solutions (with the parameters an = i (^b = 22), and the dashed curve is the 
analytical solution of this paper. The dotted vertical line indicates which of two times (T 
and T') applies. The dot-dashed line is the decay time for the Petrovay & Moreno-Insertis 
(1997) model (our equation (|26l) ). 
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Fig. 4.— Sunspot radius squared versus time, normalized by the decay time T 
numerical solution, for Bq = 7 (black) and Bq = 14 (red). Other parameters 
Figure 1. The solid curves show the numerical solution, and the dashed curves 
analytical solution. 
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Petrovay & Moreno-Insert is (1997), that the inward speed hg of the cnrrent sheet surronnding 
the decaying spot is constant. We have shown that the assnmption is jnstihed if the initial 
snnspot magnetic field Bq is not too large. Eqnation (150]) predicts that a linear decay law 
should become more accurate as Bq increases. The numerical solutions in Figure 0] confirm 
this prediction, although, as noted by the referee, they also show that the deviation from a 
linear decay is systematically underestimated in the analytical model. 

Application of the turbulent erosion theory to sunspot and starspot decay is a topic 
of current research interest (e.g., Strassmeier 2009; Rempel & Cheung 2014; Bradshaw & 
Hartigan 2014), and so our quantitative analytical predictions, reinforced by numerical so¬ 
lutions, should be useful in studies of solar and stellar activity. The value of an analytical 
calculation is that it can be used to verify more detailed magnetohydrodynamic simulations 
(e.g., Hurlburt & DeRosa 2008; Rempel & Cheung 2014) and to guide empirical models (e.g., 
Gafeira et al. 2014). 

The erosion model can be further refined. For instance, we assumed Dq = const in our 
analysis of Section 2. The rate of relative diffusion of two photospheric magnetic fragments 
is controlled by turbulent eddies whose size is equal to the current distance between the 
fragments. Consequently, the turbulent diffusivity is expected to be scale-dependent. In 
practice the turbulent diffusivity is determined by applying the induction equation to pairs 
of solar magnetograms (e.g., Chae et al. 2008, and references therein). Scale-dependent 
turbulent diffusivity has been invoked to interpret observations of photospheric flux cancel¬ 
lation (Litvinenko 2011) and the dispersion of photospheric bright points (Abramenko et al. 
2011). The turbulent erosion model of sunspot decay should be generalized to incorporate 
the dependence of the effective diffusivity on the size of a decaying sunspot. In addition, 
although Petrovay & Moreno-Insertis (1997) argued that regular radial flows play little if any 
role in sunspot decay, the effect of regular photospheric flows on sunspot decay should be 
investigated in more detail. Finally, recent observations emphasized the difference between 
the maximum and average sunspot magnetic held strengths (Tlatov & Pevtsov 2014), and 
so it may be worthwhile to derive a solution for a more general initial prohle of the magnetic 
held within the sunspot, as well as a more realistic dependence of the turbulent dihusivity 
on the held strength within the spot. 


The authors thank the referee for comments and suggestions that helped to improve the 
original manuscript. 
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A. Numerical method 


The numerical solutions in Section |3] use the Crank-Nicolson method to solve the non¬ 
linear diffusion equation (|1]), in which a discrete version of the equation is linearised at each 
time step. The Crank-Nicolson method is a preferred one for solution of parabolic partial 
differential equations because it is unconditionally stable, and second order accurate in time 
(e.g. Press et ah 1992). 


Equation (jl]) is solved at spatial locations rj = (j — l)h with j = and 

h = Tra/{L — 1), for a sequence of times = {n — l)r, with n = 1,2,.... Introducing the 
notation = B{rj,tn) and = D{B^), we consider a Crank-Nicolson scheme with the 
differencing of terms in equation (jl]): 


dB 

'~dt 


^n+l _ 


tn I'f'j 


and 


A 

dr 


rD{B)f 


tn,rj 


1 ^ 
2 dr 

1 




^n+l j 


1 

2 dr 




tn,rj 
n+1 jjn+l ' 


r ,, 1 


2h I ■’+5 


jjn+l ^ ^ 


h 


l(r. 

2h \ -^+2 i+2 h 


^-5 j-i h 

B- - BU 


2 


n+1 3 

1 ■ 

2 


h 


(Al) 


(A2) 


In the final expression in equation (IA2D . the centered differences are taken about locations 
r-i and r -, i. We introduce the approximations and 

^ 1 (C” + cy.). (A3) 

involving a linearisation in time and a spatial averaging respectively. Combining equa¬ 
tions (lAip and flA2p we have 

with s = t/K^ and 


HB';) = (j - |)£'”+syi - [(j - \)D]^ + (i - i)By] s,” + (j - DByap,. (A5) 

A von Neumann analysis of equations flA4|) - flA5D in the linear case Dj = Dq = const confirms 
that the scheme is unconditionally stable. The corresponding explicit scheme with the same 
spatial differencing is unstable if D^r/h? > ^ (e.g. Press et al. 1992). 
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Equations flA4j) and flASj) define the update for points j = 2, 3,..., L — 1. At the point 
j = 1, the boundary condition dB/dr\^^Q = 0 is enforced using the one-sided second order 
difference approximation to the derivative: 


or 


dB 

dr 


-3S" 


n+1 


+ 45?+^ 


BT' 


^n + l 


2h 


= 0 , 


(A6) 


- + 4S”+^ - 5”+^ = 0. 


(A7) 


For the point j = L we obtain an update equation allowing flux transport across the 

boundary r = via a discretisation of equation (jl]) at time t = tn and spatial location 

r = r 1 with differencing schemes 
^“2 


and 


dB 






L-i 


tn,rj 1 


r 


(AS) 


A 

dr 


rD 


.dB 

dr 


inr. 





tl-iDU 


tn,rL 


tlD" 


Bl_,-ABl_, + 3Bl 
2h 



(A9) 


where equation (lA9p involves the one-sided second order difference approximation to the 
derivative: 


dB 

dr 


TDn 

^L-2 


4:Bl^ + 3B 


tn,rL 


2h 


(AlO) 


Equations flASp and flA9D give the update equation for j = L: 


gn+l 


+ Bltl = 0 


L-2 

r _ 3 
^ 2 


D1-1 + 


3^L 


nn 

^L-2 


L-2 

1 - S^ ^ 3 ^L -1 + 3s— 


+ ( 1 - Ask — Id^ 


T — ^ 
^ 2 


\di 


Bl. 


BU 


(All) 


Equations flA4D . flASD . flA7p . and flAlip provide a system of linear equations for the held 
values with j = 1,2,..., L, which must be solved at each time step. The scheme may 

be written in matrix form as 


(h + a_A) 6*"+^ = (!' + a+A') B", 


(A12) 





























(A13) 


where B"' = (i?”, i ?2 , • • •, i a± = is the L x L matrix 

I' = diag(0,l,...,l), 

the matrix A is defined by = 0 for all j except 


All_i — —, 


a_ 


and the matrix A' is defined by A'j^j = 0 for all j except 


A' - — 

-^LL-2 ~ 


a+ \L 


L-2 


- n” -I- 
1 ^L-l + 


L - 1 


3 I 5 


^LL-1 ~ ( 1 


L - 1 


J^L I ) 


^LL = 


a+ 


. -h o , ^ u 


L - 


L - 


3 


(A14) 


(A15) 


A simple test for the new method is provided by the exact solntion for constant diffusivity 
Dq with a Ganssian profile: 

B{r, t) = ^ exp (-|r^/cr^) (A16) 

with 

cT^ = 2Z)ot + o'q, (A17) 

where the constant cxo defines the initial width. The magnetic flux (divided by 27r) from 
r = 0 to r = for this solution is 

rrm 

/ rB{r,t)dr 

Jo (A18) 

= $0 [l - exp (-Ir^/cT^)] . 

Equation flAlSp provides a check on the implementation of the boundary condition at r = r^- 
The method incurs truncation error (proportional to and at each time step, and the 
accumulation of the error limits the accuracy of the solution when the system is evolved 
over many time steps. The calculations presented in this paper are checked by trials with 
different spatial steps. 
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